Does philopatry pay? Meta-analysis reveals little empirical support for dispersal being costly.

Supplmentary Material

Author

Martinig, A. R., S. L. P. Burk, Y. Yang, M. Lagisz, and S. Nakagawa

Published

August 22, 2024

Things to do for April

  1. Rerun all the models
  2. Put figures at the end
  3. Add some texts so the reader can understand
  4. Go through code and add more annotations
  5. Ask questions about things you do not understand

1 Setting up

1.1 Load packages

Code
pacman::p_load(tidyverse, # tidy family and related pacakges below
               kableExtra, # nice tables
               MuMIn,  # multi-model inference
               emmeans, # post-hoc tests
               gridExtra, # may not use this
               pander,   # nice tables
               metafor,  # package for meta-analysis
               ape,      # phylogenetic analysis
               MuMIn,  # multi-model inference
               patchwork,   # putting ggplots together - you need to install via devtool
               here,         # making reading files easy
               orchaRd # plotting orchard plots
               # more too add

)

eval(metafor:::.MuMIn)

1.2 Load data: paper and tree data

Code
#the dataset with effect sizes
dat <- read.csv(here("data", "clean_data.csv"))

#dim(dat)
#head(dat)

#the phylogenetic tree
tree <- read.tree(here("data", "species_tree.tre"))

1.3 Calculate effect sizes and sampling variances using custom functions

Code
# function to calculate effect sizes
# Zr - correlation
# there is always n

effect_size <- function(m1, m2, sd1, sd2, n1, n2, n, # 12 arguments
                        est , se, p_val, direction, method){

  if(method == "mean_method"){
  
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s_pool <- sqrt( ((n1-1)*sd1^2 + (n2-1)*sd2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r    
    
  }else if(method == "count_method"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- sqrt(m1)
    s2 <- sqrt(m2)
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r     
    
  }else if(method == "percent_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "percent_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    m1 <- m1/100
    m2 <- m2/100
    sd1 <- sd1/100
    sd2 <- sd2/100
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method1"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(8)
    s2 <- 1/sqrt(8)
    m1 <- asin(sqrt(m1))
    m2 <- asin(sqrt(m2))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "proportion_method2"){
    
    h <- n/n1 + n/n2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    s1 <- 1/sqrt(sd1^2/(4*m1*(1-m1)))
    s2 <- 1/sqrt(sd2^2/(4*m2*(1-m2)))
    m1 <- asin(sqrt(m1/100))
    m2 <- asin(sqrt(m2/100))
    s_pool <- sqrt( ((n1-1)*s1^2 + (n2-1)*s2^2) / (n - 2) )
    j <- 1 - (3 / (4*n - 9))
    d <- ((m2 - m1) / s_pool) * j
    r_pb <-  d / sqrt(d^2 + h)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r 
    
  }else if(method == "t_method1"){
    

    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "t_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- est/sqrt(est^2 + n - 2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
    
  }else if(method == "F_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
  
  }else if(method == "F_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    r_pb <- sqrt(est)/sqrt(est + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b = r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "p_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- qt(1 - p_val, n - 2)
    r_pb <- t/sqrt(t^2 + n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r_b <- r_b*(direction)
    r <- r_b
    
  }else if(method == "correlation_method1"){
    
    r <- est
    
  }else if(method == "correlation_method2"){
    
    r <- 2*sin((pi/6)*est)
    
  }else if(method == "estimate_method1"){
    
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  }else if(method == "estimate_method2"){
    
    n1 <- n/2
    n2 <- n/2
    p <- n1/n # prop for n1
    q <- n2/n # prop for n2
    t <- est/se
    r_pb <- t/sqrt(t^2+ n -2)
    r_b <- r_pb*(sqrt(p*q)/dnorm(qnorm(p)))
    r <- r_b #r_b = r
  
  } 
  
    if(r >= 1){
    # if over 1, we use 0.99
    Zr <- atanh(0.99)

    }else if(r <= -1){

    Zr <- atanh(-0.99) # r = correlation

    } else {

    Zr <- atanh(r) # r = correlation

    }
  
  VZr <- 1 /(n - 3)
  
  data.frame(ri = r, yi = Zr , vi = VZr)
  
}


##########

#' Title: Contrast name generator
#'
#' @param name: a vector of character strings
cont_gen <- function(name) {
  combination <- combn(name, 2)
  name_dat <- t(combination)
  names <- paste(name_dat[, 1], name_dat[, 2], sep = "-")
  return(names)
}

#' @title get_pred1: intercept-less model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred1 <- function (model, mod = " ") {
  name <- firstup(as.character(stringr::str_replace(row.names(model$beta), mod, "")))
  len <- length(name)
  
   if (len != 1) {
        newdata <- diag(len)
        pred <- metafor::predict.rma(model, 
                                     newmods = newdata,
                                     tau2.levels = 1:len)
    }
    else {
        pred <- metafor::predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title get_pred2: normal model
#' @description Function to get CIs (confidence intervals) and PIs (prediction intervals) from rma objects (metafor)
#' @param model: rma.mv object 
#' @param mod: the name of a moderator 
get_pred2 <- function (model, mod = " ") {
  name <- as.factor(str_replace(row.names(model$beta), 
                                paste0("relevel", "\\(", mod,", ref = name","\\)"),""))
  len <- length(name)
  
  if(len != 1){
  newdata <- diag(len)
  pred <- predict.rma(model, intercept = FALSE, newmods = newdata[ ,-1])
  }
  else {
    pred <- predict.rma(model)
  }
  estimate <- pred$pred
  lowerCL <- pred$ci.lb
  upperCL <- pred$ci.ub 
  lowerPR <- pred$cr.lb
  upperPR <- pred$cr.ub 
  
  table <- tibble(name = factor(name, levels = name, labels = name), estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = model$pval,
                  lowerPR = lowerPR, upperPR = upperPR)
}

#' @title mr_results
#' @description Function to put results of meta-regression and its contrasts
#' @param res1: data frame 1
#' @param res1: data frame 2
mr_results <- function(res1, res2) {
  restuls <-tibble(
    `Fixed effect` = c(as.character(res1$name), cont_gen(res1$name)),
    Estimate = c(res1$estimate, res2$estimate),
    `Lower CI [0.025]` = c(res1$lowerCL, res2$lowerCL),
    `Upper CI  [0.975]` = c(res1$upperCL, res2$upperCL),
    `P value` = c(res1$pval, res2$pval),
    `Lower PI [0.025]` = c(res1$lowerPR, res2$lowerPR),
    `Upper PI  [0.975]` = c(res1$upperPR, res2$upperPR),
  )
}


#' @title all_models
#' @description Function to take all possible models and get their results
#' @param model: intercept-less model
#' @param mod: the name of a moderator 

all_models <- function(model, mod = " ", type = "normal") {
  
  # getting the level names out
  level_names <- levels(factor(model$data[[mod]]))
  dat2 <- model$data
  mod <- mod

  # creating a function to run models
  run_rma1 <- function(name) {
      
    # variance covarinace matrix for sampling errors
    VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
    VCV1 <- nearPD(VCV1)$mat
    
      
    rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }

    run_rma2 <- function(name) {
    
            # variance covarinace matrix for sampling errors
            VCVa <- vcalc(vi = dat2$vi_abs, cluster = dat$shared_group, rho = 0.5)
            VCVa <- nearPD(VCVa)$mat
               
               rma.mv(abs_yi2, V = VCVa,
               mods = ~relevel(dat2[[mod]], ref = name),
                      random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
    }
    
    run_rma3 <- function(name) {
      
              # variance covarinace matrix for sampling errors
                VCV1 <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
                VCV1 <- nearPD(VCV1)$mat
               
                rma.mv(yi, V = VCV1,
                   mods = ~relevel(dat2[[mod]], ref = name),
                     random = list( # hetero in relation to mod
                            formula(paste("~", mod, "| effectID")),
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned),
                     struct = "DIAG",
                     data = dat2,
                     test = "t",
                     sparse = TRUE,
                     control = list(optimizer = "Nelder-Mead"))
   }


# results of meta-regression including all contrast results; taking the last level out ([-length(level_names)])
# use different specifications for aboslute and hetero models 
if (type == "normal"){

    model_all <- purrr::map(level_names[-length(level_names)], run_rma1)

  }else if(type == "absolute"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma2)
    
  }else if(type == "hetero"){
    model_all <- purrr::map(level_names[-length(level_names)], run_rma3)
    
  }
  
  # getting estimates from intercept-less models (means for all the groups)
  res1 <- get_pred1(model, mod = mod)
  
  # getting estiamtes from all contrast models
  res2_pre <- purrr::map(model_all, ~ get_pred2(.x, mod = mod))
  
  # a list of the numbers to take out unnecessary contrasts
  contra_list <- Map(seq, from=1, to=1:(length(level_names) - 1))
  res2 <- purrr::map2_dfr(res2_pre, contra_list, ~.x[-(.y), ]) 
  # creating a table
  res_tab <- mr_results(res1, res2) %>% 
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")
  
  # results
  res_tab

}

1.4 Preparing final dataset

Code
# adding effect sizes to dataset
effect2 <- pmap_dfr(list(dat$mean_group_1, dat$mean_group_2, 
                         dat$variance_group_1, dat$variance_group_2, 
                         dat$n_group_1, dat$n_group_2, dat$n, 
                         dat$effect_size_value, dat$effect_size_variance, 
                         dat$effect_size_p_value_numeric, 
                         dat$direction_change, factor(dat$function_needed)), 
                    effect_size)                    

# merging two data frames
dat <- cbind(dat, effect2)

# renaming X to effectID
colnames(dat)[colnames(dat) == "X"] <- "effectID"

# creating the phylogeny column

dat$phylogeny <-  gsub(" ", "_", dat$species_cleaned)

# checking species names match between tree and dataset

# match(unique(dat$phylogeny), tree$tip.label)
# match(tree$tip.label, unique(dat$phylogeny))
# 
# intersect(unique(dat$phylogeny), tree$tip.label)
# setdiff(unique(dat$phylogeny), tree$tip.label)
# # 
# match(unique(dat$phylogeny), tree$tip.label)
# sum(is.na(match(unique(dat$phylogeny), tree$tip.label)))

# creating a correlation matrix from a tree

tree <- compute.brlen(tree)
cor_tree <- vcv(tree,corr=T) 
# creating a VCV for sampling error variances 
# we assure shared groups will have r = 0.5
# it is fine if these are not positive definite but changing this
VCV <- vcalc(vi = dat$vi, cluster = dat$shared_group, rho = 0.5)
# focing the matrix to be positive definite
# this should be in the method
VCV <- nearPD(VCV)$mat

2 Intercept meta-analytic model

2.1 Phylogenetic multilevel models

Code
# takes a while to run
mod <- rma.mv(yi = yi, 
              V = VCV,
              mod = ~ 1,
              data = dat,
              random = list(
                            ~ 1 | effectID,
                            ~ 1 | paperID,
                            ~ 1 | species_cleaned,
                            ~ 1 | phylogeny),
              R= list(phylogeny = cor_tree),
              test = "t",
              sparse = TRUE)

# saving the runs
saveRDS(mod, here("Rdata", "mod.RDS"))
Code
# getting the saved model
mod <- readRDS(here("Rdata", "mod.RDS"))

summary(mod)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.0913   558.1826   568.1826   590.7488   568.2724   

Variance Components:

            estim    sqrt  nlvls  fixed           factor    R 
sigma^2.1  0.1105  0.3324    675     no         effectID   no 
sigma^2.2  0.0110  0.1050    202     no          paperID   no 
sigma^2.3  0.0287  0.1694    146     no  species_cleaned   no 
sigma^2.4  0.0000  0.0000    146     no        phylogeny  yes 

Test for Heterogeneity:
Q(df = 674) = 390155258.0414, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0301  0.0303  -0.9947  674  0.3202  -0.0896  0.0294    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# Phylogeney accounts for nothing (0.00%), so we can take it out. 
round(i2_ml(mod),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.74              71.90               7.17              18.67 
      I2_phylogeny 
              0.00 
Code
pred_mod <- predict.rma(mod)

  estimate <- pred_mod$pred
  lowerCL <- pred_mod$ci.lb
  upperCL <- pred_mod$ci.ub 
  lowerPR <- pred_mod$cr.lb
  upperPR <- pred_mod$cr.ub 
  
  table_mod <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod <- table_mod %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_sex as RDS
saveRDS(tab_mod, here("Rdata", "tab_mod.RDS"))
Code
tab_mod <-readRDS(here("Rdata", "tab_mod.RDS"))

tab_mod
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.03 -0.09 0.029 0.32 -0.793 0.733
Code
# We remove phylogeny as a random effect from our meta-analytic model because phylogeny accounts for nothing (0%)
mod1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)

summary(mod1)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-279.0913   558.1826   566.1826   584.2355   566.2424   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1105  0.3324    675     no         effectID 
sigma^2.2  0.0110  0.1050    202     no          paperID 
sigma^2.3  0.0287  0.1694    146     no  species_cleaned 

Test for Heterogeneity:
Q(df = 674) = 390155258.0414, p-val < .0001

Model Results:

estimate      se     tval   df    pval    ci.lb   ci.ub    
 -0.0301  0.0303  -0.9947  674  0.3202  -0.0896  0.0294    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(i2_ml(mod1),2)   
          I2_Total        I2_effectID         I2_paperID I2_species_cleaned 
             97.74              71.90               7.17              18.67 
Code
orchard_plot(mod1, xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
pred_mod1 <- predict.rma(mod1)

  estimate <- pred_mod1$pred
  lowerCL <- pred_mod1$ci.lb
  upperCL <- pred_mod1$ci.ub 
  lowerPR <- pred_mod1$cr.lb
  upperPR <- pred_mod1$cr.ub 
  
  table_mod1 <- tibble("Fixed effect" = "intercept", estimate = estimate,
                  lowerCL = lowerCL, upperCL = upperCL,
                  pval = mod1$pval,
                  lowerPR = lowerPR, upperPR = upperPR)

# creating a table
tab_mod1 <- table_mod1 %>%
  kable("html",  digits = 3) %>%
  kable_styling("striped", position = "left") %>%
  scroll_box(width = "100%")

# saving tab_sex as RDS
saveRDS(tab_mod1, here("Rdata", "tab_mod1.RDS"))
Code
tab_mod1 <-readRDS(here("Rdata", "tab_mod1.RDS"))

tab_mod1
Fixed effect estimate lowerCL upperCL pval lowerPR upperPR
intercept -0.03 -0.09 0.029 0.32 -0.793 0.733

3 Uni-moderator meta-regression models

Code
#I want to move this up
#For each of our moderators, we run a uni-moderator meta-regression model
mod_class <- rma.mv(yi = yi, V = VCV,
               mod = ~ species_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_class)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-276.6580   553.3159   571.3159   611.8680   571.5891   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1090  0.3302    675     no         effectID 
sigma^2.2  0.0079  0.0891    202     no          paperID 
sigma^2.3  0.0398  0.1995    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 669) = 389349816.4547, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 669) = 0.5729, p-val = 0.7521

Model Results:

                             estimate      se     tval   df    pval    ci.lb 
species_classActinopterygii    0.0009  0.1757   0.0052  669  0.9959  -0.3441 
species_classArachnida         0.0393  0.2905   0.1353  669  0.8924  -0.5311 
species_classAves             -0.0275  0.0376  -0.7322  669  0.4643  -0.1012 
species_classInsecta           0.0943  0.1510   0.6242  669  0.5327  -0.2023 
species_classMammalia         -0.0705  0.0469  -1.5041  669  0.1330  -0.1626 
species_classReptilia          0.0672  0.1456   0.4612  669  0.6448  -0.2188 
                              ci.ub    
species_classActinopterygii  0.3459    
species_classArachnida       0.6097    
species_classAves            0.0462    
species_classInsecta         0.3908    
species_classMammalia        0.0215    
species_classReptilia        0.3531    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_class)
   R2_marginal R2_conditional 
    0.00661093     0.30922703 
Code
orchard_plot(mod_class, mod = "species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_class <- all_models(mod_class,  mod = "species_class", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_class, here("Rdata", "tab_class.RDS"))
Code
tab_class <-readRDS(here("Rdata", "tab_class.RDS"))

tab_class
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Actinopterygii 0.035 -0.240 0.309 0.805 -0.615 0.684
Arachnida 0.046 -0.391 0.484 0.835 -0.687 0.780
Aves -0.016 -0.067 0.035 0.540 -0.607 0.575
Insecta 0.135 -0.104 0.374 0.267 -0.500 0.771
Mammalia -0.014 -0.076 0.047 0.648 -0.606 0.578
Reptilia 0.147 -0.079 0.372 0.201 -0.484 0.777
Actinopterygii-Arachnida 0.048 -0.535 0.631 0.871 -0.893 0.990
Actinopterygii-Aves -0.017 -0.326 0.291 0.912 -0.818 0.784
Actinopterygii-Insecta 0.116 -0.292 0.523 0.577 -0.728 0.960
Actinopterygii-Mammalia -0.032 -0.342 0.279 0.842 -0.833 0.770
Actinopterygii-Reptilia 0.084 -0.299 0.468 0.667 -0.749 0.917
Arachnida-Aves -0.065 -0.568 0.437 0.798 -0.959 0.828
Arachnida-Insecta 0.068 -0.501 0.636 0.815 -0.865 1.000
Arachnida-Mammalia -0.080 -0.583 0.424 0.756 -0.974 0.815
Arachnida-Reptilia 0.036 -0.516 0.588 0.898 -0.886 0.958
Aves-Insecta 0.133 -0.146 0.413 0.349 -0.657 0.923
Aves-Mammalia -0.014 -0.105 0.077 0.759 -0.759 0.730
Aves-Reptilia 0.101 -0.143 0.345 0.415 -0.677 0.880
Insecta-Mammalia -0.147 -0.429 0.134 0.304 -0.938 0.643
Insecta-Reptilia -0.032 -0.393 0.329 0.862 -0.854 0.791
Mammalia-Reptilia 0.116 -0.129 0.360 0.355 -0.663 0.894
Code
mod_type <- rma.mv(yi = yi,  
               V = VCV,
               mod = ~ study_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_type)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.6859   557.3717   567.3717   589.9304   567.4617   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1103  0.3321    675     no         effectID 
sigma^2.2  0.0102  0.1008    202     no          paperID 
sigma^2.3  0.0311  0.1763    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 390147685.0977, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 673) = 0.6477, p-val = 0.5236

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
study_typenatural        -0.0283  0.0311  -0.9087  673  0.3639  -0.0894  0.0328 
study_typesemi-natural   -0.0619  0.0684  -0.9052  673  0.3657  -0.1962  0.0724 
                          
study_typenatural         
study_typesemi-natural    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_type)
   R2_marginal R2_conditional 
   0.000694167    0.272761411 
Code
orchard_plot(mod_type, mod = "study_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_type <- all_models(mod_type,  mod = "study_type", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_type, here("Rdata", "tab_type.RDS"))
Code
tab_type <-readRDS(here("Rdata", "tab_type.RDS"))

tab_type
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Natural -0.008 -0.063 0.047 0.775 -0.741 0.725
Semi-natural -0.034 -0.153 0.086 0.582 -0.774 0.707
Natural-Semi-natural -0.026 -0.143 0.092 0.671 -0.766 0.715
Code
mod_design <- rma.mv(yi = yi, V = VCV,
               mod = ~ study_design - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_design)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.8372   557.6744   567.6744   590.2332   567.7644   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1104  0.3323    675     no         effectID 
sigma^2.2  0.0118  0.1088    202     no          paperID 
sigma^2.3  0.0283  0.1684    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 390112207.7113, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 673) = 0.5268, p-val = 0.5908

Model Results:

                      estimate      se     tval   df    pval    ci.lb   ci.ub 
study_designcontrast   -0.0274  0.0319  -0.8587  673  0.3908  -0.0902  0.0353 
study_designgradient   -0.0394  0.0462  -0.8522  673  0.3944  -0.1301  0.0514 
                        
study_designcontrast    
study_designgradient    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_design)
   R2_marginal R2_conditional 
  0.0001737059   0.2670124425 
Code
orchard_plot(mod_design, mod = "study_design", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_design <- all_models(mod_design,  mod = "study_design", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_design, here("Rdata", "tab_design.RDS"))
Code
tab_design <-readRDS(here("Rdata", "tab_design.RDS"))

tab_design
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Contrast 0.007 -0.033 0.047 0.735 -0.575 0.589
Gradient -0.039 -0.106 0.028 0.257 -0.624 0.546
Contrast-Gradient -0.025 -0.107 0.057 0.553 -0.758 0.709
Code
mod_disp1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ dispersal_type - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp1)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.2694   556.5387   568.5387   595.6003   568.6650   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1107  0.3327    675     no         effectID 
sigma^2.2  0.0120  0.1095    202     no          paperID 
sigma^2.3  0.0275  0.1658    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 672) = 388370107.8889, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 672) = 0.6209, p-val = 0.6016

Model Results:

                        estimate      se     tval   df    pval    ci.lb   ci.ub 
dispersal_typeB          -0.0918  0.0757  -1.2125  672  0.2258  -0.2406  0.0569 
dispersal_typebreeding   -0.0129  0.0446  -0.2899  672  0.7720  -0.1005  0.0746 
dispersal_typenatal      -0.0273  0.0331  -0.8247  672  0.4098  -0.0924  0.0377 
                          
dispersal_typeB           
dispersal_typebreeding    
dispersal_typenatal       

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp1)
   R2_marginal R2_conditional 
   0.002031891    0.264386405 
Code
orchard_plot(mod_disp1, mod = "dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_disp1 <- all_models(mod_disp1,  mod = "dispersal_type", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_disp1, here("Rdata", "tab_disp1.RDS"))
Code
tab_disp1 <-readRDS(here("Rdata", "tab_disp1.RDS"))

tab_disp1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.042 -0.200 0.116 0.602 -0.790 0.706
Breeding 0.014 -0.065 0.094 0.727 -0.722 0.750
Natal -0.017 -0.076 0.041 0.563 -0.751 0.717
B-Breeding 0.056 -0.111 0.223 0.510 -0.694 0.806
B-Natal 0.025 -0.134 0.183 0.760 -0.724 0.773
Breeding-Natal -0.031 -0.109 0.046 0.427 -0.767 0.704
Code
mod_disp2 <- rma.mv(yi = yi, V = VCV,
               mod = ~ dispersal_phase - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_disp2)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.4741   556.9483   566.9483   589.5070   567.0382   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1103  0.3321    675     no         effectID 
sigma^2.2  0.0095  0.0973    202     no          paperID 
sigma^2.3  0.0316  0.1777    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 389177250.1681, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 673) = 1.0358, p-val = 0.3555

Model Results:

                            estimate      se     tval   df    pval    ci.lb 
dispersal_phaseprospection   -0.0666  0.0464  -1.4352  673  0.1517  -0.1576 
dispersal_phasesettlement    -0.0219  0.0320  -0.6864  673  0.4927  -0.0847 
                             ci.ub    
dispersal_phaseprospection  0.0245    
dispersal_phasesettlement   0.0408    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_disp2)
   R2_marginal R2_conditional 
   0.002008579    0.272708686 
Code
orchard_plot(mod_disp2, mod = "dispersal_phase", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_disp2 <- all_models(mod_disp2,  mod = "dispersal_phase", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_disp2, here("Rdata", "tab_disp2.RDS"))
Code
tab_disp2 <-readRDS(here("Rdata", "tab_disp2.RDS"))

tab_disp2
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Dispersal -0.047 -0.114 0.019 0.162 -0.641 0.546
Settlement 0.008 -0.034 0.051 0.698 -0.582 0.599
Dispersal-Settlement 0.048 -0.033 0.129 0.244 -0.689 0.785
Code
mod_sex <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_sex)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.2650   556.5301   568.5301   595.5916   568.6564   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1101  0.3318    675     no         effectID 
sigma^2.2  0.0118  0.1088    202     no          paperID 
sigma^2.3  0.0288  0.1698    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 672) = 389114255.4670, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 672) = 0.9049, p-val = 0.4383

Model Results:

      estimate      se     tval   df    pval    ci.lb   ci.ub    
sexB   -0.0734  0.0449  -1.6356  672  0.1024  -0.1615  0.0147    
sexF   -0.0131  0.0354  -0.3705  672  0.7112  -0.0827  0.0564    
sexM   -0.0166  0.0366  -0.4536  672  0.6503  -0.0884  0.0552    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_sex)
   R2_marginal R2_conditional 
   0.003266959    0.272185747 
Code
orchard_plot(mod_sex, mod = "sex", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_sex <- all_models(mod_sex,  mod = "sex", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_sex, here("Rdata", "tab_sex.RDS"))
Code
tab_sex <-readRDS(here("Rdata", "tab_sex.RDS"))

tab_sex
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.056 -0.139 0.027 0.184 -0.792 0.680
F 0.010 -0.053 0.074 0.749 -0.724 0.744
M 0.000 -0.066 0.066 0.994 -0.734 0.734
B-F 0.066 -0.021 0.154 0.135 -0.670 0.803
B-M 0.056 -0.034 0.145 0.222 -0.681 0.793
F-M -0.011 -0.072 0.050 0.732 -0.744 0.723

3.1 Higher life history grouping

Code
mod_age1 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age1)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-276.1033   552.2066   564.2066   591.2681   564.3329   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1105  0.3324    675     no         effectID 
sigma^2.2  0.0094  0.0970    202     no          paperID 
sigma^2.3  0.0291  0.1705    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 672) = 385949649.3241, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 672) = 2.2996, p-val = 0.0762

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
age_class_cleanadult      -0.0278  0.0315  -0.8823  672  0.3779  -0.0897 
age_class_cleanjuvenile    0.0256  0.0502   0.5093  672  0.6107  -0.0730 
age_class_cleanmix        -0.1544  0.0655  -2.3583  672  0.0186  -0.2830 
                           ci.ub    
age_class_cleanadult      0.0341    
age_class_cleanjuvenile   0.1241    
age_class_cleanmix       -0.0258  * 

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age1)
   R2_marginal R2_conditional 
   0.009541693    0.265463686 
Code
orchard_plot(mod_age1, mod = "age_class_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_age1 <- all_models(mod_age1,  mod = "age_class_clean", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_age1, here("Rdata", "tab_age1.RDS"))
Code
tab_age1 <-readRDS(here("Rdata", "tab_age1.RDS"))

tab_age1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Adult -0.004 -0.061 0.052 0.879 -0.735 0.726
Juvenile 0.012 -0.080 0.105 0.794 -0.722 0.747
Mix -0.101 -0.224 0.022 0.108 -0.839 0.638
Adult-Juvenile 0.017 -0.072 0.106 0.713 -0.717 0.751
Adult-Mix -0.096 -0.217 0.025 0.118 -0.835 0.642
Juvenile-Mix -0.113 -0.254 0.027 0.115 -0.855 0.629

3.2 Finer life stage grouping

Code
# age_class

mod_age2 <- rma.mv(yi = yi,
               V = VCV,
               mod = ~ age_class - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_age2)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-274.5087   549.0175   569.0175   614.0604   569.3523   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1109  0.3330    675     no         effectID 
sigma^2.2  0.0100  0.0999    202     no          paperID 
sigma^2.3  0.0279  0.1672    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 668) = 383749965.8761, p-val < .0001

Test of Moderators (coefficients 1:7):
F(df1 = 7, df2 = 668) = 1.3219, p-val = 0.2370

Model Results:

              estimate      se     tval   df    pval    ci.lb    ci.ub    
age_classA     -0.0189  0.0344  -0.5501  668  0.5824  -0.0865   0.0486    
age_classJ      0.0227  0.0507   0.4475  668  0.6546  -0.0768   0.1222    
age_classJA    -0.2883  0.1839  -1.5679  668  0.1174  -0.6493   0.0727    
age_classJY    -0.2084  0.0976  -2.1346  668  0.0332  -0.4000  -0.0167  * 
age_classJYA   -0.0767  0.0926  -0.8282  668  0.4079  -0.2586   0.1052    
age_classY     -0.0520  0.0642  -0.8113  668  0.4175  -0.1780   0.0739    
age_classYA    -0.0561  0.0501  -1.1193  668  0.2634  -0.1544   0.0423    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_age2)
   R2_marginal R2_conditional 
    0.01444672     0.26552449 
Code
orchard_plot(mod_age2, mod = "age_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
# table not created for this (not meaningful for some levels)

3.3 Higher level

Code
mod_fit1 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_higher_level - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit1)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.2067   556.4135   568.4135   595.4750   568.5398   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1105  0.3324    675     no         effectID 
sigma^2.2  0.0113  0.1061    202     no          paperID 
sigma^2.3  0.0285  0.1687    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 672) = 389078252.9618, p-val < .0001

Test of Moderators (coefficients 1:3):
F(df1 = 3, df2 = 672) = 0.9624, p-val = 0.4100

Model Results:

                                      estimate      se     tval   df    pval 
fitness_higher_levelindirect fitness   -0.2368  0.1706  -1.3878  672  0.1656 
fitness_higher_levelreproduction       -0.0185  0.0335  -0.5536  672  0.5800 
fitness_higher_levelsurvival           -0.0379  0.0345  -1.0971  672  0.2730 
                                        ci.lb   ci.ub    
fitness_higher_levelindirect fitness  -0.5719  0.0982    
fitness_higher_levelreproduction      -0.0842  0.0472    
fitness_higher_levelsurvival          -0.1056  0.0299    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit1)
   R2_marginal R2_conditional 
   0.002755885    0.266457669 
Code
orchard_plot(mod_fit1, mod = "fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
tab_fit1 <- all_models(mod_fit1,  mod = "fitness_higher_level", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_fit1, here("Rdata", "tab_fit1.RDS"))
Code
tab_fit1 <-readRDS(here("Rdata", "tab_fit1.RDS"))

tab_fit1
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Indirect fitness -0.226 -0.557 0.104 0.180 -1.026 0.574
Reproduction -0.005 -0.065 0.056 0.878 -0.736 0.727
Survival -0.008 -0.070 0.055 0.814 -0.739 0.724
Indirect fitness-Reproduction 0.222 -0.111 0.554 0.191 -0.579 1.022
Indirect fitness-Survival 0.219 -0.113 0.551 0.197 -0.582 1.020
Reproduction-Survival -0.003 -0.063 0.057 0.926 -0.734 0.728

3.4 Finer level

  • need to think what to do with classes less than 5 studies
Code
mod_fit2 <- rma.mv(yi = yi, 
               V = VCV,
               mod = ~ fitness_metric_clean - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_fit2)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.1855   554.3710   576.3710   625.9017   576.7741   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1105  0.3324    675     no         effectID 
sigma^2.2  0.0086  0.0925    202     no          paperID 
sigma^2.3  0.0338  0.1837    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 667) = 386811649.1407, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 667) = 0.6990, p-val = 0.6926

Model Results:

                                                           estimate      se 
fitness_metric_cleanage at first reproduction               -0.0354  0.1060 
fitness_metric_cleanindirect fitness                        -0.2312  0.1716 
fitness_metric_cleanlifetime breeding success               -0.0065  0.0391 
fitness_metric_cleanlifetime reproductive success           -0.0326  0.0426 
fitness_metric_cleanoffspring reproduction                   0.0375  0.1069 
fitness_metric_cleanoffspring survival                       0.0028  0.0481 
fitness_metric_cleanreproductive lifespan and/or attempts   -0.0682  0.1025 
fitness_metric_cleansurvival                                -0.0591  0.0385 
                                                              tval   df    pval 
fitness_metric_cleanage at first reproduction              -0.3337  667  0.7387 
fitness_metric_cleanindirect fitness                       -1.3471  667  0.1784 
fitness_metric_cleanlifetime breeding success              -0.1654  667  0.8687 
fitness_metric_cleanlifetime reproductive success          -0.7651  667  0.4445 
fitness_metric_cleanoffspring reproduction                  0.3510  667  0.7257 
fitness_metric_cleanoffspring survival                      0.0573  667  0.9543 
fitness_metric_cleanreproductive lifespan and/or attempts  -0.6651  667  0.5062 
fitness_metric_cleansurvival                               -1.5350  667  0.1253 
                                                             ci.lb   ci.ub    
fitness_metric_cleanage at first reproduction              -0.2434  0.1727    
fitness_metric_cleanindirect fitness                       -0.5681  0.1058    
fitness_metric_cleanlifetime breeding success              -0.0833  0.0704    
fitness_metric_cleanlifetime reproductive success          -0.1163  0.0511    
fitness_metric_cleanoffspring reproduction                 -0.1724  0.2475    
fitness_metric_cleanoffspring survival                     -0.0917  0.0972    
fitness_metric_cleanreproductive lifespan and/or attempts  -0.2694  0.1331    
fitness_metric_cleansurvival                               -0.1347  0.0165    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_fit2)
   R2_marginal R2_conditional 
   0.006512715    0.281583273 
Code
orchard_plot(mod_fit2, mod = "fitness_metric_clean", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
# table not created for this (not meaningful for some levels)
Code
mod_gen <- rma.mv(yi = yi, 
                V = VCV,
               mod = ~ whose_fitness - 1,
               data = dat, 
               random = list(
                 ~ 1 | effectID,
                 ~ 1 | paperID,
                 ~ 1 | species_cleaned),
               test = "t",
               sparse = TRUE)
summary(mod_gen)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.5551   557.1103   567.1103   589.6690   567.2002   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1099  0.3315    675     no         effectID 
sigma^2.2  0.0087  0.0934    202     no          paperID 
sigma^2.3  0.0338  0.1838    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 389746171.6702, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 673) = 1.1516, p-val = 0.3167

Model Results:

                         estimate      se     tval   df    pval    ci.lb 
whose_fitnessindividual   -0.0385  0.0313  -1.2292  673  0.2194  -0.0999 
whose_fitnessoffspring     0.0054  0.0464   0.1172  673  0.9067  -0.0856 
                          ci.ub    
whose_fitnessindividual  0.0230    
whose_fitnessoffspring   0.0964    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_gen)
   R2_marginal R2_conditional 
   0.001737253    0.280190472 
Code
orchard_plot(mod_gen, mod = "whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4)

Code
tab_gen <- all_models(mod_gen,  mod = "whose_fitness", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_gen, here("Rdata", "tab_gen.RDS"))
Code
tab_gen <-readRDS(here("Rdata", "tab_gen.RDS"))

tab_gen
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
Individual -0.021 -0.077 0.035 0.467 -0.756 0.714
Offspring 0.033 -0.054 0.121 0.455 -0.705 0.771
Individual-Offspring 0.054 -0.026 0.134 0.184 -0.683 0.791

3.5 Normal analysis

Code
mod_focus <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_focus)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.8197   557.6393   567.6393   590.1981   567.7293   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1103  0.3321    675     no         effectID 
sigma^2.2  0.0117  0.1080    202     no          paperID 
sigma^2.3  0.0289  0.1701    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 389708009.9183, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 673) = 0.5558, p-val = 0.5739

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0185  0.0464  -0.3988  673  0.6902  -0.1096  0.0726    
fitness_main_focusY   -0.0337  0.0319  -1.0541  673  0.2922  -0.0964  0.0291    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus)
   R2_marginal R2_conditional 
  0.0002786564   0.2693466804 
Code
orchard_plot(mod_focus, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
tab_focus <- all_models(mod_focus,  mod = "fitness_main_focus", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_focus, here("Rdata", "tab_focus.RDS"))
Code
tab_focus <-readRDS(here("Rdata", "tab_focus.RDS"))

tab_focus
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
N 0.001 -0.082 0.084 0.987 -0.734 0.736
Y -0.013 -0.069 0.044 0.655 -0.745 0.720
N-Y -0.014 -0.094 0.067 0.740 -0.748 0.721

3.6 Comparing normal model to heteroscedastic model

Code
# homo
mod_focus1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ fitness_main_focus - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_focus1)

Multivariate Meta-Analysis Model (k = 675; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-279.7147  9565.8891   569.4294   592.0029   569.5191   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1105  0.3324    675     no         effectID 
sigma^2.2  0.0115  0.1070    202     no          paperID 
sigma^2.3  0.0276  0.1662    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 389708009.9183, p-val < .0001

Test of Moderators (coefficients 1:2):
F(df1 = 2, df2 = 673) = 0.5367, p-val = 0.5849

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
fitness_main_focusN   -0.0182  0.0461  -0.3960  673  0.6922  -0.1087  0.0722    
fitness_main_focusY   -0.0328  0.0317  -1.0360  673  0.3006  -0.0950  0.0294    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_focus1)
   R2_marginal R2_conditional 
   0.000258706    0.261512359 
Code
# hetero
mod_focus2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ fitness_main_focus,
                 data = dat, 
                 random = list(
                   ~ fitness_main_focus | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_focus2)

Multivariate Meta-Analysis Model (k = 675; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-279.5317  9565.5232   571.0634   598.1517   571.1891   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0117  0.1082    202     no          paperID 
sigma^2.2  0.0266  0.1632    146     no  species_cleaned 

outer factor: effectID           (nlvls = 675)
inner factor: fitness_main_focus (nlvls = 2)

            estim    sqrt  k.lvl  fixed  level 
tau^2.1    0.1034  0.3216    162     no      N 
tau^2.2    0.1129  0.3361    513     no      Y 

Test for Residual Heterogeneity:
QE(df = 673) = 389708009.9183, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 673) = 0.1088, p-val = 0.7416

Model Results:

                     estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt               -0.0175  0.0453  -0.3872  673  0.6987  -0.1065  0.0714    
fitness_main_focusY   -0.0145  0.0439  -0.3299  673  0.7416  -0.1006  0.0717    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
#r2_ml(mod_focusA1)

# they are not sig. different
anova(mod_focus1, mod_focus2)

        df      AIC      BIC     AICc    logLik    LRT   pval             QE 
Full     6 571.0634 598.1517 571.1891 -279.5317               389708009.9183 
Reduced  5 569.4294 592.0029 569.5191 -279.7147 0.3660 0.5452 389708009.9183 
Code
# homo
orchard_plot(mod_focus1, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
# hetero
orchard_plot(mod_focus2, mod = "fitness_main_focus", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

3.7 Normal analysis

Code
mod_confirm <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_confirm)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.2069   556.4137   570.4137   601.9751   570.5827   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1106  0.3326    675     no         effectID 
sigma^2.2  0.0179  0.1337    202     no          paperID 
sigma^2.3  0.0205  0.1433    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 671) = 389396688.0027, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 671) = 0.3596, p-val = 0.8374

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0123  0.0845   0.1458  671  0.8841  -0.1536 
confirmation_biasCost      -0.0423  0.0404  -1.0459  671  0.2960  -0.1217 
confirmation_biasMixed     -0.0280  0.0400  -0.7011  671  0.4835  -0.1066 
confirmation_biasNeutral   -0.0088  0.0474  -0.1866  671  0.8520  -0.1019 
                           ci.ub    
confirmation_biasBenefit  0.1782    
confirmation_biasCost     0.0371    
confirmation_biasMixed    0.0505    
confirmation_biasNeutral  0.0842    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm)
   R2_marginal R2_conditional 
   0.001455991    0.258842550 
Code
orchard_plot(mod_confirm, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
tab_confirm <- all_models(mod_confirm,  mod = "confirmation_bias", type = "normal")

# saving tab_sex as RDS
saveRDS(tab_confirm, here("Rdata", "tab_confirm.RDS"))
Code
tab_confirm <-readRDS(here("Rdata", "tab_confirm.RDS"))

tab_confirm
Fixed effect Estimate Lower CI [0.025] Upper CI [0.975] P value Lower PI [0.025] Upper PI [0.975]
B -0.018 -0.182 0.146 0.831 -0.767 0.731
C -0.076 -0.158 0.006 0.070 -0.811 0.659
N -0.071 -0.156 0.014 0.101 -0.807 0.665
B-C -0.058 -0.229 0.113 0.504 -0.808 0.692
B-N -0.053 -0.223 0.117 0.538 -0.803 0.697
C-N 0.005 -0.093 0.103 0.923 -0.732 0.742

3.8 Comparing normal model to heteroscedastic model

Code
# homo
mod_confirm1 <- rma.mv(yi = yi, 
                V = VCV,
                mod = ~ confirmation_bias - 1,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                method = "ML",
                sparse = TRUE)
summary(mod_confirm1)

Multivariate Meta-Analysis Model (k = 675; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-279.5141  9565.4879   573.0282   604.6311   573.1961   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1110  0.3331    675     no         effectID 
sigma^2.2  0.0166  0.1288    202     no          paperID 
sigma^2.3  0.0192  0.1385    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 671) = 389396688.0027, p-val < .0001

Test of Moderators (coefficients 1:4):
F(df1 = 4, df2 = 671) = 0.3488, p-val = 0.8449

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
confirmation_biasBenefit    0.0123  0.0834   0.1475  671  0.8828  -0.1515 
confirmation_biasCost      -0.0413  0.0399  -1.0364  671  0.3004  -0.1196 
confirmation_biasMixed     -0.0264  0.0393  -0.6721  671  0.5018  -0.1035 
confirmation_biasNeutral   -0.0080  0.0467  -0.1712  671  0.8641  -0.0997 
                           ci.ub    
confirmation_biasBenefit  0.1761    
confirmation_biasCost     0.0370    
confirmation_biasMixed    0.0507    
confirmation_biasNeutral  0.0837    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_confirm1)
   R2_marginal R2_conditional 
   0.001441935    0.244972752 
Code
# hetero
mod_confirm2 <- rma.mv(yi = yi, 
                 V = VCV,
                 mod = ~ confirmation_bias,
                 data = dat, 
                 random = list(
                   ~ confirmation_bias | effectID,
                   ~ 1 | paperID,
                   ~ 1 | species_cleaned),
                 struct = "DIAG",
                 method = "ML",
                 test = "t",
                 sparse = TRUE,
                  control=list(optimizer="nlminb", rel.tol=1e-8))
summary(mod_confirm2)

Multivariate Meta-Analysis Model (k = 675; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-240.7118  9487.8833   501.4235   546.5707   501.7549   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0157  0.1251    202     no          paperID 
sigma^2.2  0.0113  0.1062    146     no  species_cleaned 

outer factor: effectID          (nlvls = 675)
inner factor: confirmation_bias (nlvls = 4)

            estim    sqrt  k.lvl  fixed    level 
tau^2.1    0.0298  0.1726     32     no  Benefit 
tau^2.2    0.1376  0.3709    204     no     Cost 
tau^2.3    0.0633  0.2516    289     no    Mixed 
tau^2.4    0.2048  0.4525    150     no  Neutral 

Test for Residual Heterogeneity:
QE(df = 671) = 389396688.0027, p-val < .0001

Test of Moderators (coefficients 2:4):
F(df1 = 3, df2 = 671) = 0.5031, p-val = 0.6802

Model Results:

                          estimate      se     tval   df    pval    ci.lb 
intrcpt                     0.0251  0.0591   0.4256  671  0.6705  -0.0908 
confirmation_biasCost      -0.0782  0.0648  -1.2077  671  0.2276  -0.2054 
confirmation_biasMixed     -0.0594  0.0614  -0.9682  671  0.3333  -0.1799 
confirmation_biasNeutral   -0.0468  0.0713  -0.6567  671  0.5116  -0.1868 
                           ci.ub    
intrcpt                   0.1411    
confirmation_biasCost     0.0490    
confirmation_biasMixed    0.0611    
confirmation_biasNeutral  0.0932    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
# they are not sig. different
anova(mod_confirm1, mod_confirm2)

        df      AIC      BIC     AICc    logLik     LRT   pval             QE 
Full    10 501.4235 546.5707 501.7549 -240.7118                389396688.0027 
Reduced  7 573.0282 604.6311 573.1961 -279.5141 77.6046 <.0001 389396688.0027 
Code
# homo
orchard_plot(mod_confirm1, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

Code
# hetero
orchard_plot(mod_confirm2, mod = "confirmation_bias", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 90)

4 Interaction meta-regression models (2 moderators)

Code
# the interaction between 2 categorical variables are the same as creating a new variable with 2 categories combined so we use that approach 

dat$sex_species_class <- paste(dat$sex, dat$species_class ,sep = "_")

mod_int1 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_species_class - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_int1)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-273.0119   546.0237   582.0237   662.8840   583.0908   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1102  0.3320    675     no         effectID 
sigma^2.2  0.0111  0.1053    202     no          paperID 
sigma^2.3  0.0315  0.1775    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 660) = 383037947.3355, p-val < .0001

Test of Moderators (coefficients 1:15):
F(df1 = 15, df2 = 660) = 0.7657, p-val = 0.7166

Model Results:

                                   estimate      se     tval   df    pval 
sex_species_classB_Aves             -0.0009  0.0584  -0.0161  660  0.9872 
sex_species_classB_Insecta          -0.1463  0.2765  -0.5291  660  0.5969 
sex_species_classB_Mammalia         -0.1797  0.0675  -2.6634  660  0.0079 
sex_species_classB_Reptilia          0.1576  0.2741   0.5750  660  0.5655 
sex_species_classF_Actinopterygii   -0.0688  0.2077  -0.3311  660  0.7407 
sex_species_classF_Arachnida         0.0402  0.2822   0.1424  660  0.8868 
sex_species_classF_Aves             -0.0361  0.0445  -0.8116  660  0.4173 
sex_species_classF_Insecta           0.2397  0.1872   1.2802  660  0.2009 
sex_species_classF_Mammalia         -0.0185  0.0538  -0.3440  660  0.7309 
sex_species_classF_Reptilia          0.0710  0.1571   0.4519  660  0.6515 
sex_species_classM_Actinopterygii    0.0683  0.2072   0.3298  660  0.7417 
sex_species_classM_Aves             -0.0346  0.0449  -0.7703  660  0.4414 
sex_species_classM_Insecta          -0.1427  0.3793  -0.3762  660  0.7069 
sex_species_classM_Mammalia         -0.0162  0.0576  -0.2817  660  0.7782 
sex_species_classM_Reptilia         -0.0241  0.2328  -0.1035  660  0.9176 
                                     ci.lb    ci.ub     
sex_species_classB_Aves            -0.1157   0.1138     
sex_species_classB_Insecta         -0.6891   0.3966     
sex_species_classB_Mammalia        -0.3121  -0.0472  ** 
sex_species_classB_Reptilia        -0.3806   0.6959     
sex_species_classF_Actinopterygii  -0.4767   0.3391     
sex_species_classF_Arachnida       -0.5139   0.5942     
sex_species_classF_Aves            -0.1235   0.0513     
sex_species_classF_Insecta         -0.1279   0.6073     
sex_species_classF_Mammalia        -0.1242   0.0872     
sex_species_classF_Reptilia        -0.2375   0.3796     
sex_species_classM_Actinopterygii  -0.3385   0.4752     
sex_species_classM_Aves            -0.1226   0.0535     
sex_species_classM_Insecta         -0.8876   0.6022     
sex_species_classM_Mammalia        -0.1294   0.0969     
sex_species_classM_Reptilia        -0.4811   0.4330     

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_int1)
   R2_marginal R2_conditional 
    0.01827497     0.29189297 
Code
orchard_plot(mod_int1, mod = "sex_species_class", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
dat$sex_whose_fitness <- paste(dat$sex, dat$whose_fitness ,sep = "_")

mod_int2 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_whose_fitness - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_int2)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.0899   554.1798   572.1798   612.7318   572.4529   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1102  0.3319    675     no         effectID 
sigma^2.2  0.0088  0.0938    202     no          paperID 
sigma^2.3  0.0332  0.1821    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 669) = 386931618.7208, p-val < .0001

Test of Moderators (coefficients 1:6):
F(df1 = 6, df2 = 669) = 0.9616, p-val = 0.4504

Model Results:

                               estimate      se     tval   df    pval    ci.lb 
sex_whose_fitnessB_individual   -0.0784  0.0456  -1.7181  669  0.0862  -0.1679 
sex_whose_fitnessB_offspring     0.0200  0.1506   0.1330  669  0.8943  -0.2756 
sex_whose_fitnessF_individual   -0.0114  0.0372  -0.3057  669  0.7599  -0.0844 
sex_whose_fitnessF_offspring    -0.0209  0.0550  -0.3805  669  0.7037  -0.1288 
sex_whose_fitnessM_individual   -0.0338  0.0382  -0.8851  669  0.3764  -0.1087 
sex_whose_fitnessM_offspring     0.0612  0.0631   0.9694  669  0.3327  -0.0628 
                                ci.ub    
sex_whose_fitnessB_individual  0.0112  . 
sex_whose_fitnessB_offspring   0.3157    
sex_whose_fitnessF_individual  0.0617    
sex_whose_fitnessF_offspring   0.0870    
sex_whose_fitnessM_individual  0.0412    
sex_whose_fitnessM_offspring   0.1852    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_int2)
   R2_marginal R2_conditional 
   0.006784346    0.280770163 
Code
orchard_plot(mod_int2, mod = "sex_whose_fitness", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
dat$sex_fitness_higher_level <- paste(dat$sex, dat$fitness_higher_level ,sep = "_")

mod_int3 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_fitness_higher_level - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_int3)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.1147   554.2295   576.2295   625.7602   576.6326   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1104  0.3323    675     no         effectID 
sigma^2.2  0.0120  0.1097    202     no          paperID 
sigma^2.3  0.0288  0.1698    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 667) = 387165630.2892, p-val < .0001

Test of Moderators (coefficients 1:8):
F(df1 = 8, df2 = 667) = 0.7222, p-val = 0.6720

Model Results:

                                            estimate      se     tval   df 
sex_fitness_higher_levelB_reproduction       -0.0230  0.0719  -0.3201  667 
sex_fitness_higher_levelB_survival           -0.0948  0.0515  -1.8388  667 
sex_fitness_higher_levelF_indirect fitness   -0.2311  0.2549  -0.9064  667 
sex_fitness_higher_levelF_reproduction       -0.0038  0.0391  -0.0977  667 
sex_fitness_higher_levelF_survival           -0.0258  0.0440  -0.5862  667 
sex_fitness_higher_levelM_indirect fitness   -0.2299  0.2079  -1.1059  667 
sex_fitness_higher_levelM_reproduction       -0.0280  0.0428  -0.6539  667 
sex_fitness_higher_levelM_survival           -0.0007  0.0452  -0.0149  667 
                                              pval    ci.lb   ci.ub    
sex_fitness_higher_levelB_reproduction      0.7490  -0.1642  0.1182    
sex_fitness_higher_levelB_survival          0.0664  -0.1960  0.0064  . 
sex_fitness_higher_levelF_indirect fitness  0.3651  -0.7316  0.2695    
sex_fitness_higher_levelF_reproduction      0.9222  -0.0806  0.0730    
sex_fitness_higher_levelF_survival          0.5580  -0.1121  0.0606    
sex_fitness_higher_levelM_indirect fitness  0.2692  -0.6381  0.1783    
sex_fitness_higher_levelM_reproduction      0.5134  -0.1121  0.0561    
sex_fitness_higher_levelM_survival          0.9881  -0.0895  0.0881    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_int3)
   R2_marginal R2_conditional 
   0.007263485    0.275448328 
Code
orchard_plot(mod_int3, mod = "sex_fitness_higher_level", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, angle = 45)

Code
dat$sex_dispersal_type <- paste(dat$sex, dat$dispersal_type ,sep = "_")

mod_int4 <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ sex_dispersal_type - 1,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              test = "t",
              sparse = TRUE)
summary(mod_int4)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-277.2697   554.5393   578.5393   632.5548   579.0171   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1106  0.3325    675     no         effectID 
sigma^2.2  0.0156  0.1247    202     no          paperID 
sigma^2.3  0.0248  0.1576    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 666) = 387281534.0754, p-val < .0001

Test of Moderators (coefficients 1:9):
F(df1 = 9, df2 = 666) = 0.4575, p-val = 0.9029

Model Results:

                              estimate      se     tval   df    pval    ci.lb 
sex_dispersal_typeB_B          -0.0961  0.1147  -0.8371  666  0.4029  -0.3214 
sex_dispersal_typeB_breeding   -0.0041  0.0997  -0.0409  666  0.9674  -0.1998 
sex_dispersal_typeB_natal      -0.0834  0.0511  -1.6323  666  0.1031  -0.1838 
sex_dispersal_typeF_B          -0.0812  0.1119  -0.7255  666  0.4684  -0.3008 
sex_dispersal_typeF_breeding   -0.0041  0.0539  -0.0755  666  0.9398  -0.1099 
sex_dispersal_typeF_natal      -0.0080  0.0395  -0.2014  666  0.8405  -0.0855 
sex_dispersal_typeM_B          -0.1084  0.1232  -0.8802  666  0.3790  -0.3504 
sex_dispersal_typeM_breeding   -0.0217  0.0587  -0.3693  666  0.7120  -0.1370 
sex_dispersal_typeM_natal      -0.0041  0.0415  -0.0981  666  0.9219  -0.0855 
                               ci.ub    
sex_dispersal_typeB_B         0.1293    
sex_dispersal_typeB_breeding  0.1917    
sex_dispersal_typeB_natal     0.0169    
sex_dispersal_typeF_B         0.1385    
sex_dispersal_typeF_breeding  0.1017    
sex_dispersal_typeF_natal     0.0696    
sex_dispersal_typeM_B         0.1335    
sex_dispersal_typeM_breeding  0.0936    
sex_dispersal_typeM_natal     0.0774    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
r2_ml(mod_int4)
   R2_marginal R2_conditional 
   0.006429907    0.272231265 
Code
orchard_plot(mod_int4, mod = "sex_dispersal_type", xlab = "Effect Size: Zr", group = "paperID", branch.size = 4, trunk.size = 0.5, angle = 45)

5 Multi-moderator-regression model

Code
# use "ML" for model selection
#######################
# Mulit-variable models
#######################

# we need - to
mod_full <- rma.mv(yi = yi, 
               V = VCV,
              mod = ~ 1 + sex +
                age_class_clean + whose_fitness + 
                fitness_higher_level +
                dispersal_type + dispersal_phase,
              data = dat, 
              random = list(
                ~ 1 | effectID,
                ~ 1 | paperID,
                ~ 1 | species_cleaned),
              method = "ML",
              test = "t",
              sparse = TRUE)
summary(mod_full)

Multivariate Meta-Analysis Model (k = 675; method: ML)

   logLik   Deviance        AIC        BIC       AICc   
-273.8369  9554.1335   575.6737   638.8797   576.3101   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1088  0.3299    675     no         effectID 
sigma^2.2  0.0082  0.0906    202     no          paperID 
sigma^2.3  0.0309  0.1757    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 664) = 381605233.2172, p-val < .0001

Test of Moderators (coefficients 2:11):
F(df1 = 10, df2 = 664) = 1.2033, p-val = 0.2853

Model Results:

                                  estimate      se     tval   df    pval 
intrcpt                            -0.3090  0.1869  -1.6531  664  0.0988 
sexF                                0.0436  0.0502   0.8694  664  0.3850 
sexM                                0.0421  0.0509   0.8280  664  0.4080 
age_class_cleanjuvenile             0.1054  0.0562   1.8758  664  0.0611 
age_class_cleanmix                 -0.0845  0.0684  -1.2355  664  0.2171 
whose_fitnessoffspring              0.0548  0.0486   1.1268  664  0.2602 
fitness_higher_levelreproduction    0.2101  0.1705   1.2320  664  0.2184 
fitness_higher_levelsurvival        0.1784  0.1732   1.0301  664  0.3033 
dispersal_typebreeding              0.0234  0.0848   0.2763  664  0.7824 
dispersal_typenatal                 0.0043  0.0774   0.0559  664  0.9554 
dispersal_phasesettlement           0.0360  0.0536   0.6713  664  0.5023 
                                    ci.lb   ci.ub    
intrcpt                           -0.6760  0.0580  . 
sexF                              -0.0549  0.1421    
sexM                              -0.0578  0.1420    
age_class_cleanjuvenile           -0.0049  0.2158  . 
age_class_cleanmix                -0.2187  0.0498    
whose_fitnessoffspring            -0.0407  0.1502    
fitness_higher_levelreproduction  -0.1247  0.5448    
fitness_higher_levelsurvival      -0.1616  0.5184    
dispersal_typebreeding            -0.1430  0.1899    
dispersal_typenatal               -0.1477  0.1564    
dispersal_phasesettlement         -0.0693  0.1413    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_full)*100, 2)
   R2_marginal R2_conditional 
          2.01          27.91 
Code
# multi-model selection
candidates <- dredge(mod_full, trace = 2)

# saving tab_sex as RDS
saveRDS(candidates, here("Rdata", "candidates.RDS"))
Code
candidates <-readRDS(here("Rdata", "candidates.RDS"))

candidates
Global model call: rma.mv(yi = yi, V = VCV, mods = ~1 + sex + age_class_clean + 
    whose_fitness + fitness_higher_level + dispersal_type + dispersal_phase, 
    random = list(~1 | effectID, ~1 | paperID, ~1 | species_cleaned), 
    data = dat, method = "ML", test = "t", sparse = TRUE)
---
Model selection table 
   (Int) age_cls_cln dsp_phs dsp_typ ftn_hgh_lvl sex whs_ftn df   logLik  AICc
1      +                                                      4 -264.109 536.3
33     +                                                   +  5 -263.311 536.7
3      +                   +                                  5 -263.495 537.1
2      +           +                                          6 -262.717 537.6
35     +                   +                               +  6 -262.796 537.7
34     +           +                                       +  7 -262.007 538.2
17     +                                           +          6 -263.045 538.2
4      +           +       +                                  7 -262.071 538.3
9      +                                       +              6 -263.248 538.6
41     +                                       +           +  7 -262.356 538.9
49     +                                           +       +  7 -262.438 539.1
36     +           +       +                               +  8 -261.430 539.1
11     +                   +                   +              7 -262.513 539.2
5      +                           +                          6 -263.726 539.6
37     +                           +                       +  7 -262.794 539.8
18     +           +                               +          8 -261.792 539.8
19     +                   +                       +          7 -262.821 539.8
10     +           +                           +              8 -261.821 539.9
43     +                   +                   +           +  8 -261.988 540.2
25     +                                       +   +          8 -262.084 540.4
42     +           +                           +           +  9 -261.057 540.4
12     +           +       +                   +              9 -261.067 540.4
51     +                   +                       +       +  8 -262.240 540.7
7      +                   +       +                          7 -263.282 540.7
50     +           +                               +       +  9 -261.232 540.7
6      +           +               +                          8 -262.446 541.1
39     +                   +       +                       +  8 -262.456 541.1
20     +           +       +                       +          9 -261.439 541.2
38     +           +               +                       +  9 -261.543 541.4
57     +                                       +   +       +  9 -261.582 541.4
44     +           +       +                   +           + 10 -260.597 541.5
27     +                   +                   +   +          9 -261.740 541.8
21     +                           +               +          8 -262.792 541.8
26     +           +                           +   +         10 -260.822 542.0
8      +           +       +       +                          9 -261.876 542.0
13     +                           +           +              8 -262.947 542.1
52     +           +       +                       +       + 10 -260.907 542.2
45     +                           +           +           +  9 -261.981 542.2
53     +                           +               +       +  9 -262.077 542.4
40     +           +       +       +                       + 10 -261.051 542.4
15     +                   +       +           +              9 -262.341 543.0
58     +           +                           +   +       + 11 -260.343 543.1
28     +           +       +                   +   +         11 -260.363 543.1
59     +                   +                   +   +       + 10 -261.401 543.1
22     +           +               +               +         10 -261.537 543.4
14     +           +               +           +             10 -261.591 543.5
23     +                   +       +               +          9 -262.646 543.6
46     +           +               +           +           + 11 -260.649 543.7
47     +                   +       +           +           + 10 -261.719 543.8
54     +           +               +               +       + 11 -260.816 544.0
29     +                           +           +   +         10 -261.866 544.1
16     +           +       +       +           +             11 -260.880 544.2
55     +                   +       +               +       + 10 -261.961 544.3
60     +           +       +                   +   +       + 12 -260.047 544.6
24     +           +       +       +               +         11 -261.238 544.9
48     +           +       +       +           +           + 12 -260.249 545.0
61     +                           +           +   +       + 11 -261.302 545.0
56     +           +       +       +               +       + 12 -260.543 545.6
31     +                   +       +           +   +         11 -261.594 545.6
30     +           +               +           +   +         12 -260.584 545.7
62     +           +               +           +   +       + 13 -259.966 546.5
63     +                   +       +           +   +       + 12 -261.181 546.9
32     +           +       +       +           +   +         13 -260.161 546.9
64     +           +       +       +           +   +       + 14 -259.716 548.1
   delta weight
1   0.00  0.110
33  0.43  0.089
3   0.80  0.074
2   1.28  0.058
35  1.44  0.054
34  1.91  0.042
17  1.94  0.042
4   2.03  0.040
9   2.35  0.034
41  2.61  0.030
49  2.77  0.028
36  2.80  0.027
11  2.92  0.026
5   3.30  0.021
37  3.48  0.019
18  3.53  0.019
19  3.54  0.019
10  3.59  0.018
43  3.92  0.016
25  4.11  0.014
42  4.12  0.014
12  4.13  0.014
51  4.42  0.012
7   4.46  0.012
50  4.47  0.012
6   4.84  0.010
39  4.86  0.010
20  4.88  0.010
38  5.09  0.009
57  5.17  0.008
44  5.26  0.008
27  5.48  0.007
21  5.53  0.007
26  5.71  0.006
8   5.75  0.006
13  5.84  0.006
52  5.88  0.006
45  5.96  0.006
53  6.15  0.005
40  6.16  0.005
15  6.68  0.004
58  6.82  0.004
28  6.86  0.004
59  6.86  0.004
22  7.14  0.003
14  7.25  0.003
23  7.29  0.003
46  7.43  0.003
47  7.50  0.003
54  7.77  0.002
29  7.79  0.002
16  7.89  0.002
55  7.98  0.002
60  8.30  0.002
24  8.61  0.001
48  8.71  0.001
61  8.74  0.001
56  9.30  0.001
31  9.32  0.001
30  9.38  0.001
62 10.22  0.001
63 10.57  0.001
32 10.61  0.001
64 11.81  0.000
Models ranked by AICc(x) 
Code
# displays delta AICc <2
candidates_aic2 <- subset(candidates, delta < 2) 
# model averaging
mr_averaged_aic2 <- summary(model.avg(candidates, delta < 2)) 

mr_averaged_aic2

Call:
model.avg(object = candidates, subset = delta < 2)

Component model call: 
rma.mv(yi = yi, V = VCV, mods = ~<7 unique rhs>, random = list(~1 | 
     effectID, ~1 | paperID, ~1 | species_cleaned), data = dat, method = ML, 
     test = t, sparse = TRUE)

Component models: 
       df  logLik   AICc delta weight
(Null)  4 -264.11 536.28  0.00   0.24
4       5 -263.31 536.71  0.43   0.19
2       5 -263.49 537.08  0.80   0.16
1       6 -262.72 537.57  1.28   0.12
24      6 -262.80 537.72  1.44   0.11
14      7 -262.01 538.19  1.91   0.09
3       6 -263.05 538.22  1.94   0.09

Term codes: 
age_class_clean dispersal_phase             sex   whose_fitness 
              1               2               3               4 

Model-averaged coefficients:  
(full average) 
                           Estimate Std. Error z value Pr(>|z|)
intrcpt                   -0.026596   0.039875   0.667    0.505
whose_fitnessoffspring     0.020470   0.036076   0.567    0.570
dispersal_phasesettlement  0.012379   0.029609   0.418    0.676
age_class_cleanjuvenile    0.004109   0.022391   0.183    0.854
age_class_cleanmix        -0.020044   0.047793   0.419    0.675
sexF                       0.005767   0.022722   0.254    0.800
sexM                       0.005266   0.021666   0.243    0.808
 
(conditional average) 
                          Estimate Std. Error z value Pr(>|z|)
intrcpt                   -0.02660    0.03988   0.667    0.505
whose_fitnessoffspring     0.05191    0.04085   1.271    0.204
dispersal_phasesettlement  0.04553    0.04141   1.099    0.272
age_class_cleanjuvenile    0.01917    0.04528   0.423    0.672
age_class_cleanmix        -0.09351    0.06154   1.520    0.129
sexF                       0.06468    0.04450   1.453    0.146
sexM                       0.05906    0.04569   1.292    0.196
Code
# relative importance of each predictor for all the models
importance <- sw(candidates)

importance
                     whose_fitness dispersal_phase age_class_clean
Sum of weights:      0.42          0.37            0.33           
N containing models:   32            32              32           
                     fitness_higher_level sex  dispersal_type
Sum of weights:      0.24                 0.23 0.15          
N containing models:   32                   32   32          

6 Publication bias

Code
# write conditional statement here

# if each group n is avaiable - assume n/2
dat$effectN <- ifelse(is.na(dat$n_group_1), (dat$n/2)*2/dat$n,  
                      (dat$n_group_1 * dat$n_group_2) / (dat$n_group_1 + dat$n_group_2))
  
dat$sqeffectN <- sqrt(dat$effectN)

mod_egger <- rma.mv(yi = yi, 
                    V = VCV,
                mod = ~ sqeffectN,
                data = dat, 
                random = list(
                  ~ 1 | effectID,
                  ~ 1 | paperID,
                  ~ 1 | species_cleaned),
                test = "t",
                sparse = TRUE)
summary(mod_egger)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.8842   557.7684   567.7684   590.3271   567.8583   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1106  0.3326    675     no         effectID 
sigma^2.2  0.0118  0.1084    202     no          paperID 
sigma^2.3  0.0279  0.1671    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 390060408.7356, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 673) = 0.0049, p-val = 0.9443

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -0.0317  0.0404  -0.7846  673  0.4330  -0.1111  0.0477    
sqeffectN    0.0003  0.0045   0.0698  673  0.9443  -0.0085  0.0092    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_egger)*100, 2)
   R2_marginal R2_conditional 
          0.00          26.41 
Code
small <- bubble_plot(mod_egger,
                     mod = "sqeffectN",
                     group = "paperID",
                     xlab = "sqrt(Effective N)",
                     ylab = "Effect Size: Zr",
                     g = TRUE)

small

Code
# creating publication year from "reference"
dat$reference <- as.character(dat$reference)
dat$year <- with(dat, substr(reference, nchar(reference)-4, nchar(reference)))
dat$year <- as.integer(ifelse(dat$year == "2017a" | dat$year == "2017b", 2017, dat$year))
# decline effect
mod_dec <- rma.mv(yi = yi, V = vi,
                     mod = ~ year,
                     data = dat, 
                     random = list(
                       ~ 1 | effectID,
                       ~ 1 | paperID,
                       ~ 1 | species_cleaned),
                     test = "t",
                     sparse = TRUE)
summary(mod_dec)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-224.2401   448.4802   458.4802   481.0390   458.5702   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.0637  0.2524    675     no         effectID 
sigma^2.2  0.0121  0.1102    202     no          paperID 
sigma^2.3  0.0277  0.1663    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 673) = 7836.0459, p-val < .0001

Test of Moderators (coefficient 2):
F(df1 = 1, df2 = 673) = 0.5799, p-val = 0.4466

Model Results:

         estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt   -2.3721  3.0862  -0.7686  673  0.4424  -8.4318  3.6876    
year       0.0012  0.0015   0.7615  673  0.4466  -0.0018  0.0042    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_dec)*100, 2)
   R2_marginal R2_conditional 
          0.20          38.58 
Code
decline <- bubble_plot(mod_dec,
                       mod = "year",
                       group = "paperID",
                       xlab = "Publication year",
                       ylab = "Effect Size: Zr",
                       g = TRUE)
decline 

Code
# All together

mod_comb <- rma.mv(yi = yi, 
                   V = VCV,
                   mod = ~ year + sqeffectN,
                   data = dat, 
                   random = list(
                     ~ 1 | effectID,
                     ~ 1 | paperID,
                     ~ 1 | species_cleaned),
                   # ~ 1 | phylogeny),
                   # R= list(phylogeny = cor_tree),
                   test = "t",
                   sparse = TRUE)
                   
summary(mod_comb)

Multivariate Meta-Analysis Model (k = 675; method: REML)

   logLik   Deviance        AIC        BIC       AICc   
-278.4433   556.8866   568.8866   595.9481   569.0129   

Variance Components:

            estim    sqrt  nlvls  fixed           factor 
sigma^2.1  0.1102  0.3319    675     no         effectID 
sigma^2.2  0.0110  0.1048    202     no          paperID 
sigma^2.3  0.0307  0.1752    146     no  species_cleaned 

Test for Residual Heterogeneity:
QE(df = 672) = 389760666.0556, p-val < .0001

Test of Moderators (coefficients 2:3):
F(df1 = 2, df2 = 672) = 0.1674, p-val = 0.8459

Model Results:

           estimate      se     tval   df    pval    ci.lb   ci.ub    
intrcpt     -1.9785  3.3810  -0.5852  672  0.5586  -8.6172  4.6601    
year         0.0010  0.0017   0.5755  672  0.5651  -0.0023  0.0043    
sqeffectN   -0.0002  0.0046  -0.0359  672  0.9713  -0.0091  0.0088    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Code
round(r2_ml(mod_comb)*100, 2)
   R2_marginal R2_conditional 
          0.09          27.50 
Code
# small-study
bubble_plot(mod_egger,
            mod = "sqeffectN",
            group = "paperID",
            xlab = "sqrt(Effective N)",
            ylab = "Effect Size: Zr",
            g = TRUE)

Code
# decline
bubble_plot(mod_comb,
            mod = "year",
            group = "paperID",
            xlab = "Publication year",
            ylab = "Effect Size: Zr",
            g = TRUE)

Code
dat <- dat %>%
  mutate(leave_out = reference)

dat$leave_out <- as.factor(dat$leave_out)


LeaveOneOut_effectsize <- list()
for (i in 1:length(levels(dat$leave_out))) {
  temp_dat <- dat %>%
    filter(leave_out != levels(dat$leave_out)[i])
  
  VCV_leaveout <- vcalc(vi = temp_dat$vi, cluster = temp_dat$shared_group, rho = 0.5)
  
  LeaveOneOut_effectsize[[i]] <-  rma.mv(yi = yi,
                                         V = VCV_leaveout, 
                                         random = list(
                                           ~ 1 | effectID,
                                           ~ 1 | paperID,
                                           ~ 1 | species_cleaned),
                                         test = "t",
                                         sparse = TRUE,
                                         data = temp_dat[temp_dat$leave_out != levels(temp_dat$leave_out)[i], ])
}

# writing function for extracting est, ci.lb, and ci.ub from all models
est.func <- function(model) {
  df <- data.frame(est = model$b, lower = model$ci.lb, upper = model$ci.ub)
  return(df)
}

# using dplyr to form data frame
MA_oneout <- lapply(LeaveOneOut_effectsize,function(x) est.func(x)) %>%
  bind_rows %>%
  mutate(left_out = levels(dat$leave_out))


# telling ggplot to stop reordering factors
MA_oneout$left_out <- as.factor(MA_oneout$left_out)
MA_oneout$left_out <- factor(MA_oneout$left_out, levels = MA_oneout$left_out)

# saving the runs
saveRDS(MA_oneout, here("Rdata", "MA_oneout.RDS"))
Code
MA_oneout <- readRDS(here("Rdata", "MA_oneout.RDS"))

# plotting
leaveoneout <- ggplot(MA_oneout) + geom_hline(yintercept = 0, lty = 2, lwd = 1) +

  geom_hline(yintercept = mod1$ci.lb, lty = 3, lwd = 0.75, colour = "black") +
  geom_hline(yintercept = mod1$b, lty = 1, lwd = 0.75, colour = "black") + 
  geom_hline(yintercept = mod1$ci.ub,
             lty = 3, lwd = 0.75, colour = "black") + 
  geom_pointrange(aes(x = left_out, y = est,
                      ymin = lower, ymax = upper)) + 
  xlab("Study left out") + 
  ylab("Zr (effect size), 95% CI") +
  coord_flip() + 
  theme(panel.grid.minor = element_blank()) + theme_bw() + theme(panel.grid.major = element_blank()) +
  theme(panel.grid.minor.x = element_blank()) + theme(axis.text.y = element_text(size = 6))

leaveoneout

8 R Session Information

Code
# pander for making it look nicer
sessionInfo() %>% pander()

R version 4.2.0 (2022-04-22)

Platform: x86_64-apple-darwin17.0 (64-bit)

locale: en_US.UTF-8||en_US.UTF-8||en_US.UTF-8||C||en_US.UTF-8||en_US.UTF-8

attached base packages: stats, graphics, grDevices, utils, datasets, methods and base

other attached packages: orchaRd(v.2.0), here(v.1.0.1), patchwork(v.1.1.3), ape(v.5.7-1), metafor(v.4.4-0), numDeriv(v.2016.8-1.1), metadat(v.1.2-0), Matrix(v.1.5-4.1), pander(v.0.6.5), gridExtra(v.2.3), emmeans(v.1.8.9), MuMIn(v.1.47.5), kableExtra(v.1.3.4), lubridate(v.1.9.3), forcats(v.1.0.0), stringr(v.1.5.1), dplyr(v.1.1.3), purrr(v.1.0.2), readr(v.2.1.4), tidyr(v.1.3.0), tibble(v.3.2.1), ggplot2(v.3.4.4) and tidyverse(v.2.0.0)

loaded via a namespace (and not attached): httr(v.1.4.7), jsonlite(v.1.8.7), viridisLite(v.0.4.2), splines(v.4.2.0), highr(v.0.10), stats4(v.4.2.0), vipor(v.0.4.7), yaml(v.2.3.7), pillar(v.1.9.0), lattice(v.0.22-5), glue(v.1.6.2), digest(v.0.6.33), rvest(v.1.0.3), colorspace(v.2.1-0), sandwich(v.3.1-0), htmltools(v.0.5.6.1), pkgconfig(v.2.0.3), xtable(v.1.8-4), mvtnorm(v.1.2-3), scales(v.1.3.0), webshot(v.0.5.5), svglite(v.2.1.2), tzdb(v.0.4.0), timechange(v.0.2.0), mgcv(v.1.9-0), farver(v.2.1.1), generics(v.0.1.3), TH.data(v.1.1-2), pacman(v.0.5.1), withr(v.2.5.2), cli(v.3.6.1), survival(v.3.5-7), magrittr(v.2.0.3), estimability(v.1.4.1), evaluate(v.0.23), fansi(v.1.0.5), nlme(v.3.1-163), MASS(v.7.3-60), xml2(v.1.3.5), beeswarm(v.0.4.0), tools(v.4.2.0), hms(v.1.1.3), lifecycle(v.1.0.4), multcomp(v.1.4-25), munsell(v.0.5.0), compiler(v.4.2.0), systemfonts(v.1.0.5), rlang(v.1.1.1), grid(v.4.2.0), rstudioapi(v.0.15.0), htmlwidgets(v.1.6.1), labeling(v.0.4.3), rmarkdown(v.2.25), gtable(v.0.3.4), codetools(v.0.2-19), R6(v.2.5.1), zoo(v.1.8-12), knitr(v.1.45), fastmap(v.1.1.1), utf8(v.1.2.3), rprojroot(v.2.0.4), mathjaxr(v.1.6-0), latex2exp(v.0.9.6), ggbeeswarm(v.0.7.2), stringi(v.1.7.12), parallel(v.4.2.0), Rcpp(v.1.0.11), vctrs(v.0.6.4), tidyselect(v.1.2.0), xfun(v.0.40) and coda(v.0.19-4)